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Abstract 

Wc propose an implicit iterative algorithm for an exact penalty method arising from in- 
equality constrained optimization problems. A rapidly convergent fixed point method is 
developed for a regularized penalty functional. The applicability and feasibility of the 
proposed method is demonstrated using large scale inequality constrained problems. 
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1 Introduction 

Let us consider the constrained optimization 

min F(x), (1) 

subject to the unilateral constraint 

{Gx - g), < 0. (2) 

The bilateral constraint £i < {Gx — g)i < Ui can be transformed into the unilateral constraint, and we only 
consider ((2} without loss of generality. We assume F is a smooth functional on R" and G G R'"^^" ig onto. 

Inequality constrained optimization problems appear in a vast range of applications such as contact prob- 
lems [23], obstacle problems [12) . topology optimization [T] 1151 I26| . robotics and gait analysis [4], contact 
mechanics [27], and there are several numerical methods which can be used as practical tools for solving the 
problems; [131 1161 1171 1181 1211 1241 1281 129) . For references to the literature on the numerical methods for opti- 
mization problems, one may also refer to the monographs [6l 1111 [121 122) . 

Interior or exterior penalty methods require solving a sequence of unconstrained problems in which the 
penalty parameter (the controlling parameter) approaches or infinity. This yields the ill-conditioning of the 
unconstrained problem, which is the main drawback of the penalty method. In contrast, exact penalty methods 
transform the constrained problem H]) - ([2]) into a single unconstrained problem. Surprisingly, the penalized 
unconstrained problems are exact under certain sufficient conditions for a local optimality in the problem ((TJ - 
([2]), i.e., all solutions of the penalized unconstrained problem are also solutions of the original problem for all 
values of the penalty parameter grater than some positive value. For this reason, considerable attention has 
been devoted to the use of exact penalty approaches in solving constrained optimization problems. A survey of 
the chronological development of the penalty methods (including multiplier methods) since 1968 to 1993 can be 
found in [7]. 

The exact penalty methods require, however, minimization of a nondifferentiable cost functional. One 
may not employ a standard optimization solver that are customized for optimization problems with smooth 
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functions. Therefore numerical techniques should be developed by utilizing the particular structures of the 
penalty functions that compensate for the absence of differentiability. Numerical methods to approximate the 
solution of exact penalty methods have been considered by several authors. We only mention the articles of 
[H El [HI El 1101 1111 1141 1251 129j . In this article we consider the exact penalty formulation where max is used for 
the penalty; 

min J(x) = F{x) + /3ip(Gx-g) 
with = ^max(0,yi)> 2/ £ R™, 

i 

for P > 0, and we develop fast iterative methods for finding the minimizer based on the nonsmooth optimization 
theory. 

The optimality condition of ((Sjl is given by 

-F'{x)& pG'di,iGx-g), (4) 
where St/j is the convex sub-differential of ip, i.e., 

dil^iy) = {s G iJ™ : Si G amax(0,yO. 1 < « < m}, 



9max(0, s) — 



[0,1], s = 0, 
1, s>0. 



On the other hand, the necessary optimality of fT])-© is given by 

i^'(x-) + GV = 0, fi = ma.x{0, fi + [Gx - g)), (5) 

where jj, G R™ is the Lagrange multiplier of the unilateral constraint [22]. Let the pair {x,jl) be a solution to 
([5]). Then x is also a solution to ([4]) provided that /3 > maxi (e.g., [5]). 

Due to the singularity and the non-uniquness of the subgradient of 9^, the direct treatment of the condition 
Q many not be feasible for numerical computation. The common strategy to alleviate the technical difficulty 
resulting from the non-differentiability of the penalty functional is to introduce a regularized penalty functional: 
Let us consider the regularized problem to ([H]); 

min J^{x) = F{x) + f3-tf},{Gx ~ g) 
with = J/SR™' 

i 

where (/f>e(s) for s G R is a regularization of the function s — >■ max(0, s) defined by 

' ^ s < 



s e 
2^ 2 



s G [0, £] 



s > e. 



(7) 



An arbitrary e > in is used to avoid the singularity and to determine a single value in the subdifferential 
9max(0,s). Since (f)^ £ G^ , the necessary optimality condition of ([6]) is given by the equation 

F'{x)+pG'i'',{Gx-g)^0. (8) 

Although the non-uniqueness for concerning subdifferential in the optimality condition ((4} is now bypassed 
through regularization, the optimality condition ([SJ is still nonlinear. One of the strategies for solving ((S)) is to 
use the asymptotic solution at infinity of the nonlinear ODE 

^^-GF'ixit))-PG'i^UGx{t)^g), 

with some positive definite matrix G which serves as a precondition of F'{x{t)). (See [29]). The method is simple 
and easy to implement, however, the convergence speed is quite slow and the numerical solution obtained by the 
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algorithm is not accurate. This is due to the fact that the nonUnearity in ^'^{Gx{t) — g) is not fully taken into 
account and not incorporated in algorithms. One of the objective of this paper is to design the fast, accurate 
numerical algorithm for ([8]) by taking the nonlinearity into consideration. 

The outline of the paper is as follows. In Section 2 an implicit iterative algorithm for (jS)) is proposed. The 
property and convergence of the proposed algorithm are analyzed. In Section 3 the Primal-Dual Active method 
is introduced and the relation to the proposed method is discussed. In Section 4 several numerical tests are 
reported to assess the performance of the method. 

2 Algorithm and Convergence Analysis 

In this section we introduce the algorithm for ^ and analyze its convergence. First, we have the consistency 
result as e — >■ 0^. 

Theorem 1. Let be a solution of the regularized problem (|6]). For an arbitrary /3 > 0, any cluster point of 
{a;e}e>o IS also a solution of (pljl . 

Proof. Let x* be a solution to ([Sjl. Then we have 



As a consequence of Theorem [TJ we have that 

Corollary 1. Suppose that each of the problem ((3]) and ((6| admits the unique solution. If 13 > Q is sufficiently 
large, x^ converges to the solution x of 

2.1 Successive iteration algorithm 

We propose the fast algorithm that provides an accurate numerical solution of For this objective, we first 
observe that the necessary optimality condition ^ is written as 



F{x,) +l3-4^,{x,) < F{x*) + P^p,{x*). 



(9) 



Let a; be a cluster point of {x^}tyo- From @, we have 



F{x) + I3^{x) + F{x,) - F{x) + - ij{x,) + ij{x,) - V'(S)) 

< Fix") +I3^{x'') + PiM^l - ^(a:*)). 

Thus, from continuity of F, ^ and the fact that < tl^eix) — tp{x) < — for all x, we obtain 



F{x) + PiIj{x) < F{x*) + /3 



□ 



F'(x) + l3G\x^i^)Gx - f,{x)) = 0, 



(10) 



where Xe(^) denotes a diagonal matrix with the entries 




{Gx 
{Gx 



9) J > 



(11) 



and f^{x) is a column vector depending on x defined by 




{Gx 
{Gx 



g)j > 
g)j<o 



(12) 



The optimality condition in the form (|10|l suggests the following fixed point iteration; 

aP{x''+^ - x'') + F'{x'') + /3G'(x.(2;'=)Ga:'=+' - f,{x'')) = 0, 



(13) 
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where P is positive, symmetric and serves a pre-conditioner for F"{xk)- The parameter a > serves a stabilizing 
and acceleration stepsize (see, Theorem 3). 

Let d*" = x'''^^ — x'^. Equation JTSj for x'^^^ is equivalent to the equation for d*" 

aPd" + F'ix") + pG'{x^{x'')G{x'' + d'=+') - Mx")) = 0, 

which gives us 

{aP + PG\4x'')G)d'' = -JUx''). (14) 
Lemma 1 The direction d* is a descent direction for Jc{x) at x''. 

Proof. From (fTI)) 

(dfe, j;(a;'=)) = -((QP + /JG'x.(a;''')G)dfc,d'=) = -a(Pd\d'=) -/3(x.(2:')Gd^Gd'=) < 0, 
where we used the fact that P is strictly positive definite. □ 

So the iteration (|13|) can be seen as a descent method and is written as; 
Algorithm 1: Fixed point iteration ([13}. 

Step 0. Set parameters: /3, a, e, P. 

Step 1. Compute the direction by (P + jSG^ x^{x'')G)d'' = -J'{x''). 
Step 2. Update = x'' + d*. 

If I J'(a::'°)loo < TOL, then stop. Otherwise repeat Step 1 - Step 2. 
Let us make some remarks on the Algorithm: 

Remark 1. In many applications, the structure of F'{x), P and G are sparse block diagonals, and the resulting 
system (|14p for the direction d* then becomes a linear system with a sparse symmetric positive- definite matrix, 
and can be efficiently solved by, for example the Cholesky decomposition method. 

Remark 2. If F{x) = ^{x, Ax) — {b,x), then we have F' {x) = Ax — b. For this case we may use the alternative 
update 

aP{x''+^ - x'') + Ax''+^ -b + p G\x4x'')Gx''+^ - f,{x'')) = 0, (15) 
assuming that it doesn't cost much to perform 

iaP + A + pG'x,{x'')G)d'' = -JUx''). 

Algorithm 1 is globally and rapidly convergent and the following results justify the fact. Let us introduce 
some notations; A'' = {j : {Gx'' - g)j > 0} and l'' = {i : [{Gx'' - g)i < 0}. 

Lemma 1. Let R{x,x) :— —{F{x) — F{x) — F' {x){x — x)) + a {P{x — x),x — x), and Xk = Xt{^'^)- The following 
identity holds for all k; 

Rix''+\x'') + F{x''+') - Fix") + I (xfcGd^ Gd") + i ^ {[XkkjAiGx'^' - g),f - {{Gx"- - g),\') = 0. 

Proof Muhiplying ^ by d'' = x''+^ - x*' 

a [Pd'', d'') - {F{x''+^) ~ Fix'') - F'ix'')d'') + F(a;'''+') - Fix'') + Ek = 0, 

where 

E,=PixtiGx''+' -g),Gd''). 
From the identity 2(a, a - b) = \\a - 6|p + ||a||^ - ||6|p with a = ^Xfc(Gx''+^ - g) and b = ^^iGx'' - g), we 
obtain 

£fc = (a, a - 6) = I iXkG'd", Gd") + i ^ ilxkkj, i(Gx'=+^ - g)A' - liGx"' - g),f). 

□ 
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Theorem 2. Assume there exists uj > such that 

R(x,x) > uj\\x - x\\'^, x,xGR". 
If there exists ko such that X'^^^ C for all k > ko, then 

Lj - x'^f + J,(x'=+i) - J,(x'') < 

and {x^} is globally convergent. 

Proof. Since s — >■ ^^(-ys) on s > is concave (see (O for the definition of (j>{s)), we have 

MVt) - MVt) < {M^s)y{s -t) = — ^L^-Lis - 1), vs,t> 0, 

max(e, yt) lyt 



tlius 



Hence 



\ Y.^^^''^o,oMGx''+' - g),\-' -\(Gx'^ - g),\'')> ^ [U{Gx''+' ~ g),) - UiGx"^ - 9),)] 
= MGx''+' - 5) - MGx" -9)^Y. [^^HG^"^' - 9).) - UiGx" - ,g)0] 



>MGx''+^ -g)-MGx'' ^g). 
Thus, we obtain 

J.(a;*=+') + R{x''+\x'') + ^ (xfcG(a:'+' - x""), G(a:*^+' - x"")) < Mx"). 

If we assume R{x,x) > i^\\x — i|p for some lj > 0, then Ji{x^) is monotonically decreasing and 

^ <oo. 

fc> fcg 

□ 

Corollary 2. Suppose i G I* 6ut (Ga;''+^ - g); > 0. A ssume 

I (XfcG(a;'--+^ - a;"), Gix'^' - x")) - p MiGx^^' - g).) > -co' \\x''+' - x'f 

with < uj' < UJ, then the algorithm is globally convergent. 

Algorithm 1 closely resembles to the semismooth Newton's method [22] applied to the equation (|4]): The 
gradient ip'^iGx — g) aX x = x^ has a Newton derivative G^ Nk{Gx^ — g) where the diagonal matrix is defined 
by 

r i, (G^'=-,g),G (0,e) 
[NkU = <^ ^ (16) 
I. 0, otherwise. 

Replacing Step 1 with the system for the semi-smooth Newton step 

(F"(a;'=) + /?G*iVfeG)d!v = -J'.{x''), (17) 

one arrives at a semismooth Newton's method. In general, the sequence {a;*'} generated by the Newton's 
method is guaranteed to converge when the initial guess x'^ is sufficiently close to the true solution x*. When 
x^ is not close enough to the minimum, taking the full Newton step d% need not decrease the objective function 
Je{x), moreover it may generate a non-convergence sequence. On the other hand, through several numerical 
experiment the sequence generated by Algorithm 1 converges to the true solution within a few iterations even 
when an initial guess is far from the true solution. 
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Lemma 2. // the fc**" iterate is close to the solution x* and satisfies Gx'' 
Algorithm 1 enjoys the superUnear convergence of semi-smooth Newton's method. 



g < e, then Xk = and 



We shall investigate Algorithm 1 and the semi-smooth Newton through a simple problem: we consider the 
optimization problem 



mmmize - a- 

2" " 



- (3tp^{x - g). 



In this case (|13p and (|17p are explicitly given as 



k+l 



${x), ${x)^x- 



I max ( 3: - 

X + 7 )— 

' max (a;- 



3,0) 



1 + 7 



Bign(3:-3) 
max(3: — g,e) 



sign(s) = 



1, s > 0, 
0, s < 0. 



and 



^ I max(3:-3,0) 
' max(x — 



1 + 7 



3ign(3: — g) 



(18) 



(19) 



(20) 



It is easy to prove that the sequence {x'^} generated by the iteration (|19|l and (|20|) converges to x — g-\- {—g)j^ 
for any initial x" provided that <; < and j3 > e — g. On the other hand, if we assume that g > Q, then for any 
initial a;" we have a;*^ = for all k > 2. 

We depict y — ${x), y = $jv(x) in Fig. [1] and [2] respectively. The outcomes after three iterations starting 
from a::o = —1.2 are also plotted to visualize the iteration process. The parameters g = —1, e — 10~^ and 
P = 2 were used to draw these graphs. (We select a large e — 10~^ just for the purpose of the visualization. In 
practice, we will take much smaller number, say, e — 10~^.) We observe from the figure that the performance 
of Algorithm 1 is much better than the one of the semismooth Newton's method. 
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Figure 1: Visualization of convergence pattern 
of the successive iteration algorithm 1, Xfc+i = 



Figure 2: Visualization of convergence pattern 
of the nonsmooth Newton's method, x^+i = 



2.2 Successive iteration with the line search 

Since d*^ determined by (|14p is a descent direction of Jt, one can use the line search method; 
Algorithm 2: Successive iteration with line search 

Step 0. Set parameters : 7, e, P. 

Step 1. Compute the direction by 



(P + 7G'x'G)d'= 



-J\x^). 



(21) 
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Step 2. Determine the steplength otk by minimizing J{x^ + ad!'), i.e., a* — argminc J(x''' + ad^). 
Step 3. Update = a;'' + a'^d'''. 

If I J'(a;'-')| < TOL, then stop. Otherwise repeat Step 1 - Step 3. 

Step 2 can be replaced by the line search algorithms such as Armijo's rule [22| . The proof of the convergence of 
Algorithm 2 is quite standard and we omit the proof. 

3 Semi-smooth method 

An alternative to our gradient-based algorithm is the Newton update. The semi-smooth Newton method in [22] 
is based on complementarity condition for F' -\- G^n = 0; 

^ = max(0, ^ + 7 (Ga; — g)), V7 > 0. 

The semi-smooth Newton method reduces to the Primal-Dual Active set method [22| ; 
Primal-Dual Active set metliod 

• Choose {x^\ //') and set fc = 0. 

• Set A*" = {j : + ■y{Gx'' - g))j > 0} and = {i : K/' + 7(6*2;'^ - g)), < 0}. 

• Solve for (a:'^+\/+i) 

F'(a;'=+i)-f GV'''+^ =0, ^'=+^=0on2'= 

(Ga;'=+i - ,g)j = on ^^ 

• Convergent or set k = k + 1 and Return to Step 2. 

Remark 3. (1) If F is quadratic, i.e., F' — Ax — b, then Step 3 is written as 

A {G^y \ f x''+' \ _ f b 
G^ )[ m'^^^ )'[g 

If A > 0, then Step 3 is solvable. Otherwise, we assume that A is positive on N{{G{j, :)). 

(2) If F IS C^, the Newton step for Step 3 is given by 

- x'') + + G*A'=+i = 

(3) In general one may use the regularized update 

A {G^Y \ ( x''+'' \ _ f b 
-el ) [ M^+i )-{g 

to avoid the possible singularity of linear system H22|) . It reduces to 

Ax''+' + P {G^y^ (g''x''+' - ff) : 



(22) 



which is very similar to psp with a = 0. Consequently, Algorithm 1 is much stabler than Prima-Dual Active 
method. 

(4) But, it IS shown in ]20/ if the Primal-Dual Active method converges, it converges q-super linearly and in a 
finite step. 

(5) One can hybrid Algorithm 1 and Prima-Dual Active method so that one may accelerate the convergence. 
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4 Numerical tests 



In this section we show some numerical experiments using Algorithm 1 proposed in Section 2 for unilateral 
constrained quadratic optimization problems 

F{x)^^{x,Ax)-{x,b), Gx<g 

with several A, b, G and g. All tests confirm the fact convergence and effectiveness of the proposed algorithm. 

4.1 Example 1: Obstacle problem 

Let Q = [0, 1] X [0, 1]. We solve an obstacle problem 



- / |Vm|^ dx — C I u{x) dx 

^ Jn Jn (23) 

K = {ue Ho{n) I u{x) < S(x,dn)}, 



mm 



where C = 10 is used and 5{x,di}) is distance from x to dQ; 

5{x,dQ) = i (1 - max(l2x- - 1|, \2y - 1|)) . 

We use the standard bilinear finite element method to discretize the problem: For Cartesian grid {xi,yj) = 
(n'n)' — — '^1 *6 define a finite element by {K,Qh,M); the element domain K is a. rectangle, K = 
[xi, Xi+i] X [j/j, yj+i], and the space of shape functions Qh = Uh € L^(Q) is given by 

"^i^Ik = [1' — Ax~' Ax ' ' ' "'+^'^+^ ' '"''J+ij • 

and A/" is nodal variables at the grid points. The subscript h indicates the mesh size h = i.The finite element 
discretization yields the discrete energy functional 

^ / IVithl^ dx - C [ Uh dx = \{Uh, HhUh)m^ - [Fh, Uh)m" 
Jn Jn ^ 

for Uh £ Qh, where Uh = (^to,o,uo,i, • • • ,Mn,n) and Hu and Fh denote the stiffness matrix and the load vector 
associated with the discretization. The inequality constrained is approximated by Uh{zi,j) < S{zi,j,dQ,), Zij = 
{xi,yj), which is equivalent to Uh < g, where g = {S{zofi,d^l), . . . , 5{zn,n, dfl), —5{zofi, . . . , — 
Let {Ufi}k denote the generated sequence by Algorithm 1. We report 

MUh) = ^{utHhU^)^„2 - {fh,uj;)^^„2 +PMUh-gh). 

and the sup norm of the gradient of J^iU^): 

||VJ.((7^)||oo = max| (HhUt - h + fiVMUk - 5;.)) J- 

We run Algorithm 1 with the following parameters and preconditioner: 

1. mesh size h = 0.02, = 0.01, P = Hh, a = 1, e = . 

2. mesh size h = 0.01, P = 0.01, P = Hh, a = 1, e = . 

Fig. |3]shows the monotone convergence of the objective function J^: the convergence achieves after 11 iteration 
for h = 0.02, and 20 iteration for h = 0.01. 
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Figure 3: The monotone decreasing of J{U^) and || V J(?7|) ||oo of the finite element solutions 
{U^} generated by Algorithm 1 with /? = 0.01, P = H^, a = 1 and e = h?. 

4.2 Example 2: Inverse source identification problem 

Let Vl = [0, 1]^. The problem consists in recovering the source term u G L^(Sl) in the equation 

- Ay = ii in y = 0, on T, (24) 

from the noisy data ys{x) £ L^ip.) such that ys = y{x) + 5{x) where 5{x) is an additive (unknown) noise. We 
assume that the source term u is constrained; < u{x) < 1 for a.e. x £ Q. Let y{u) denote the (weak) solution 
of (|24p . The problem is well-known to be ill-posed and an approximation to the solution u can be obtained by 
Tikhonov regularization method: 



min F{u) = i / (y{u) - ysYdx + 



n 



u G L^(Q), < it(x) < 1 for a.e. x G S7, 



where ys G L^{Q), t; > is a regularization parameter. Algorithm 1 requires the computation of the gradient 
F'{u). One can calculate the gradient by 

F'{u) ^p + rj u, 

where the adjoint variable p is obtained by solving the adjoint equation 

-Ap = y{u) - y\ 

In our computation, the noisy data ys is generated by adding a random noise to the observation data y: 

ys{x) = y{x) +Ta.nd{x), 

where rand(a::) is a uniformly distributed random function in [—1,1], and S is the noise level. The unknown 
source (exact solution) u and the noise free data y are depicted in the first row of Fig. |4l 

The domain = [0, 1]^ is divided into 60^ subsquares of the mesh size h — The central finite difference 
method is used to approximate —Au at Xij; 



And we approximate F{u) as follows: 



^ j^iyi'u-) - ysfdx + J^u^dx ^ Y^(b(ii)]i,j - [yshjf + ^X^"' 



2 



(25) 



(26) 



Hence the discretized exact penalty problem is equivalent to 



mmJ-\\K-'u-ys\\l„2+^\\u\\l„2 + p M^)- (27) 
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Here K is the matrix for the second order central difTerence associated to ((25}. We employed Algorithm 1 to 
the problem (|27|l with parameters a = 1, P = 1, e — . The preconditioner P = K~^K~^ + rjl is used: Step 1 
in Algorithm 1 is written as 

i[K-'YK-' + /3x.(w'))d' = -iP + Pu + 

which is equivalently written as 

{I + r^K'' + l3K\,{u''))dt' = -A'2(p + ,^u'=+M("''')), 

where we use = K. The reconstructed source obtained by the nonsmooth Tikhonov regularization with 
the regularization parameter rj, and the noisy data with noise level 5 are shown in Fig. 2] We observed that 
Algorithm 1 converged (i.e., \ Je{u'')\ao < 10~^*) within 20 iterations in all cases. The regularization parameter 
r) was selected manually according to the noise level 5. The study of automated selection of r) can be found in 
vast literature on Tikhonov regularization. 

4.3 Inverse medium problem 

Consider the inverse medium problem; determine the potential function u{x) > in 

^/^y + uy^f, yeHl,{n) (28) 

from measurement ys{x) = y{x) + S{x) of the potential y. The problem can be casted as a constrained least 
square problem; find u 

min -|y(u) - zs\l2^n^ + ||u|i2(n) 

subject to < It < [/ with a priori upper bound U, where y{u) is the solution to (|28|) . One can calculate F'{u) 
using the adjoint equation 

— Ap up — y{u) — zs 

i.e., 

F (u) — —y{u)p + rju. 

In the computation, we use the function in FigO (top left) as the unknown potential to be recovered. All the 
noisy data was generated by adding a random noise to the exact data y: 

ys{x) = y{x) + (5max{j/(a::)}rand(a::), 

where y was computed by solving the equation y = (—A + u)^^/ with / = 10. The noise free data y is depicted 
in Fig[S](top right). 

As the preconditioner in Algorithm 1, we used P = + r). Since G — I, the matrix P + Pxei^^) in Step 1 
is diagonal. Hence, the computation of the decent direction df^ is cheap but one faces the slow convergence of 
the algorithm due to the poorly chosen preconditioner. More than 1000 time iteration was required to meet the 
stopping criterion |VJe(ti )joo < 10^''' in each test. 



5 Application to nonsmooth Tikhonov regularization 

The implicit fixed point iteration proposed in Section 2 is also applicable to the optimization problem involving 
4>[s) := |s| (or the sum J]]^ 4>[xi)) in tl^e objective function, for instance, 

2 2 
^ n n 

min J{u) ■.= -\\u- vec{f)\\l^. + 771 ^ (0([O.u]O + 4>{[Dyu],)) + 772 ^ 4>{[Hu],), 
"eK" ^ ,=1 ,=1 

2 

where / G R"^" is a given noisy image, vec{f) = /2,i, . . . , fn,n] G R" , Dx, Dy represent finite differences 
in X and y direction and H denotes a discrete Laplacian. The problem is obtained by discretizing the multi- 
parameter nonsmooth Tikhonov regularization for a denoising problem 

1 2 

min -\\u- /||^2 +r]x\\u\\TV + r?2 1| AuH^i . 
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Source, u 



Data, y 





Figure 4: Example 2 (Inverse source identification). Reconstructed source (left) and the noisy 
data (right). Mesh 60 x 60. h = Parameters used in Algorithm 1 are; a = 1, P = K*K + 7], 
(3 = 1 and e = h?. The regularization parameter rj was selected manually according to the 
noise level 5. 



11 



Potential, u Data, y 




Figure 5: Example 3 (Potential identification problem). Reconstructed potential (left) and 
the noisy data (right). Mesh 60 x 60. h = Parameters used in Algorithm 1 are; q = 1, 
P = + r/, /3 = 1 and e = h^. The regularization parameter r/ was selected manually 
according to the noise level 6. 
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Here r/i , 7^2 are regularization parameters which must be appropriately selected in order to obtain a desired 
reconstructed image |19l I20| . We define (p^, the regularization of <j), by 

s s > e 



Ms) = ' 



s e - - 



s e 



se [-e,0] 
s < -e. 



The derivative is written as (j)i{s) = 
the successive iteration algorithm: 



max(e, |s|) 



, and one follows the similar argument in Section 2 to arrive at 



Update: u'^^ — + d*. 



Here the diagonal matrix Xt{'") G 



for V eM." is defined by 
1 



and Je is a regularization of J; 

2 2 

n n 

j=i 1=1 
Another example that the proposed method can handle includes the denosing problem by Total variation: 



J sjul+ul 



dxdy 



Let (j>[s) be a function defined for s > by 



Ms) 



< s < e 



2e 2 

The regularized objective functional takes the form 



Ms) 



< s < e 



e < s 



>^e(w) = ^11" ~ +0^^ 4>e{\Jul +ul)dxdy 



Let -i/Je (u) be a discretization of the second term 



where = \/[Dxu\l + [Dyit]^. The derivative of t/jj («) 



dtpeiu) _sr^,i, ,[D^u\iDx{i,k) + [Dyu]iDy{i,k) 



duk 



Thus we have 



V i[D,U],D4i,k) + [DyU],Dy{i, k)) 

^ max(e, Ti) rt 
ip'^{u) = Dix^{u)D:,u + Dlxe{u)DyU 
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Here the diagonal matrix Xe(^) is defined by 



From the observation we arrive at the successive iteration algorithm for solving the nonlinear equation J'^ (u) = 0: 

(P + Dix.{u'')D^ + Dlx.{u)Dy)d'' = -U^"), =1"+ d\ 

The details of the method and the numerical tests will be reported elsewhere. 
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